A review of statistical models used to characterize species-habitat associations with movement data: Appendix 3
Katie Florko, Ron Togunov, Marie Auger-Méthé
24/01/2024
1 Overview
This document contains a coding tutorial that demonstrates how to perform the analyses associated with the case study in “A review of statistical models used to characterize species-habitat associations with animal movement data”. All analyses link the movement data of a ringed seal in eastern Hudson Bay, Canada, to modeled prey diversity.
2 Set-up
2.1 Load packages
# data wrangling
library(here)
library(tidyverse)
library(lubridate)
# mapping
library(rnaturalearth)
library(rnaturalearthdata)
library(rgdal)
#if you need to install rgdal, install from archive:
#devtools::install_version('rgdal', version = '1.6-7') )
library(terra)
library(tidyterra)
library(sf)
library(viridis)
# fitting models
library(amt) # for selection functions
library(momentuHMM) # for hidden Markov model
library(geosphere)
# plotting
library(ggplot2)2.2 Prey diversity data
Next we load in the fish data, which contains the spatial distribution of prey diversity in 2012. The data is in a .csv, and we will rasterize this using it’s specific coordinate reference system. This is a subset of the data from Florko et al. (2021a, 2021b).
# load prey diversity data
fish <- read.csv(here("data/fish.csv"))
# rasterize prey diversity data
fish_raster <- terra::rast(fish, crs = "+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")Visualize the fish data.
# prepare land data for mapping
natearth <- ne_countries(scale = "medium",returnclass = "sf")
natearth <- natearth[which(natearth$continent!="Antarctica"),]
# projection: project land to match fish raster
nat_trans <- st_transform(natearth, crs(fish_raster))
# plot fish map
fishmap <- ggplot() +
tidyterra::geom_spatraster(data = fish_raster) +
scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
geom_sf(data = nat_trans, fill = "grey80", color = "white") +
coord_sf(xlim = ext(fish_raster)[1:2], ylim = ext(fish_raster)[3:4], expand = F)
fishmap
2.3 Seal data
Next we load in one movement track from a ringed seal equipped with an ARGOS satellite telemetry transmitter over the course of over four months in the winter of 2012-2013. This is a subset of the data from Florko et al. (2023a).
# load seal data
seal <- read.csv(here("data/seal_track_m.csv"))
head(seal)## lon lat date id
## 1 357940.9 -368949.6 2012-10-29 14:00 116484_1
## 2 371594.3 -353092.9 2012-10-30 14:00 116484_1
## 3 388429.2 -361478.4 2012-10-31 14:00 116484_1
## 4 399326.2 -357288.1 2012-11-01 14:00 116484_1
## 5 411788.2 -349500.7 2012-11-02 14:00 116484_1
## 6 407708.1 -372564.6 2012-11-03 14:00 116484_1# ensure the data is in the correct format
seal <- seal %>%
mutate(id = as.character(id),
date = as.Date(date)) Visualize the seal data on top of the fish data.
# plot seal and fish data together
sealfishmap <- fishmap +
geom_point(data=seal, aes(x=lon, y=lat), alpha = 0.6, color = "#FCEEAE") +
geom_path(data=seal, aes(x=lon, y=lat), color = "#FCEEAE")
sealfishmap
3 Fit models
We will fit four main models: 1) a resource selection function (RSF), 2) a step selection function (SSF) without habitat covariates in the movement kernel, 3) a SSF with a habitat covariate modifying the movement kernel, and 4) a hidden Markov model (HMM). All four of these main models will include prey diversity as a covariate. Note that we also fit additional HMMs to demonstrate their flexibility.
All of the step selection functions (both the RSF and two SSFs) will be fit using the amt package (Signer et al. 2019), while the HMMs will use the functions from the momentuHMM package (McClintok and Michelot 2018).
3.1 RSF
The RSF is the simplest of the four models. Fitting an RSF to movement data first requires us to generate a sample of availability points and extract covariates for the used and available locations.
# prep data and generate availability sample
set.seed(2023)
data_rsf <- seal %>%
make_track(lon, lat, date) %>% # convert data to track format
random_points() # generate availability sample; default is ten times as many available points as observed points
# plot used vs available locations on-top of prey diversity
data_rsf_map <- sealfishmap +
geom_point(data=data_rsf, aes(x=x_, y=y_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type")
data_rsf_map
We can see that the availability sample is generated within the minimum convex polygon of the used samples.
We now extract covariate (prey diversity) values for the used and available locations.
data_rsf <- data_rsf %>%
extract_covariates(fish_raster)Next, we will fit the model. The response, case_, is the column in data_rsf that identifies whether the location is a used (TRUE) or available (FALSE) location, and preydiv is the covariate for prey diversity.
# fit rsf (a binomial logistic regression)
rsf1 <- data_rsf %>%
amt::fit_rsf(case_ ~ preydiv, model = TRUE)View the summary.
# see model summary
summary(rsf1)##
## Call:
## stats::glm(formula = formula, family = stats::binomial(link = "logit"),
## data = data, model = TRUE)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5753 -0.4729 -0.4294 -0.3754 2.3363
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.118 1.400 -4.369 1.25e-05 ***
## preydiv 6.353 2.312 2.748 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 938.28 on 1539 degrees of freedom
## Residual deviance: 930.60 on 1538 degrees of freedom
## AIC: 934.6
##
## Number of Fisher Scoring iterations: 5We see that prey diversity is a significant positive covariate. We do not interpret the intercept, since it is not ecologically meaningful in a RSF. See Fieberg et al. (2021) for a detailed discussion on how to interpret parameters.
3.2 SSF
Next we will fit our two SSFs. The workflow is similar to that of the RSF, however, the availability sample is generated differently. We transform the seal locations into a track format, then convert the track data into step format (i.e., with a start and an end), and then generate the availability sample.
# prep data and generate availability sample
data_obs_steps <- seal %>%
make_track(lon, lat, date) %>% # convert data to track format
steps()
set.seed(2023)
data_ssf <- data_obs_steps %>% # convert track data to step format (i.e., with a start and an end)
random_steps() # generate availability sample By default, when applied to objects of class step_xy, the function random_steps() samples the step lengths from a gamma distribution with parameters estimated from the observed step lengths and turning angles from a von Mises distribution with parameters estimated from the observed turning angles.
Visualize the sample created.
# plot used vs available locations on-top of prey diversity
data_ssf_map <- sealfishmap +
geom_point(data=data_ssf, aes(x=x2_, y=y2_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type")
data_ssf_map
We can see that the availability sample is generated at each step and is not restricted to the minimum convex polygon.
We now extract prey diversity values for the used and available locations.
# extract prey diversity covariate
data_ssf <- data_ssf %>%
extract_covariates(fish_raster, where = "end") #sample at end of stepWe will transform both the step length and turning angle using log and cosine transformations, respectively.
# transform movement covariates
data_ssf <- data_ssf %>%
mutate(cos_ta = cos(ta_),
log_sl = log(sl_))Next, we will fit the models. The response, case_, is the column of data_ssf that identifies whether the location is a used (TRUE) or available (FALSE) location, and preydiv is the covariate for prey diversity at that location. sl_ is step length, log_sl is the log transformation of step length, and cos_ta is the cosine transformation of turning angle. strata(step_id_) specifies that the this is a conditional logistic regression that groups data by step identification number.
We are fitting two SSFs: one without a movement-related covariate (called ssf1), and one with a movement-related covariate (called ssf2).
# fit ssfs
## ssf1: ssf without covariate affecting movement kernel
ssf1 <- data_ssf %>%
fit_clogit(case_ ~ preydiv + sl_ + log_sl + cos_ta + strata(step_id_), model = TRUE)
## ssf2: ssf with covariate affecting movement kernel
ssf2 <- data_ssf %>%
fit_clogit(case_ ~ sl_ + preydiv*log_sl + cos_ta + strata(step_id_), model = TRUE)First we will interpret ssf1.
summary(ssf1)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv + sl_ +
## log_sl + cos_ta + strata(step_id_), data = data, model = TRUE,
## method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## preydiv 9.475e-01 2.579e+00 6.837e+00 0.139 0.890
## sl_ -6.155e-06 1.000e+00 1.606e-05 -0.383 0.701
## log_sl 4.219e-02 1.043e+00 1.570e-01 0.269 0.788
## cos_ta 3.066e-02 1.031e+00 1.307e-01 0.235 0.814
##
## exp(coef) exp(-coef) lower .95 upper .95
## preydiv 2.579 0.3877 3.908e-06 1.702e+06
## sl_ 1.000 1.0000 1.000e+00 1.000e+00
## log_sl 1.043 0.9587 7.668e-01 1.419e+00
## cos_ta 1.031 0.9698 7.981e-01 1.332e+00
##
## Concordance= 0.516 (se = 0.03 )
## Likelihood ratio test= 0.24 on 4 df, p=1
## Wald test = 0.24 on 4 df, p=1
## Score (logrank) test = 0.24 on 4 df, p=1This model estimates the coefficients (coef) as well as the exponent of the coefficient (exp(coef)). The coefficients are as in our regular logistic regression (the RSF), and the exp(coef) quantifies the relative intensity of use of two locations that differ by 1 unit of prey diversity but are otherwise the same. The model suggests our seal would be 2.6 times more likely to choose a location with 1 unit higher prey diversity (see Fieberg et al. 2021). However, this is not significant (p = 0.890), and we see that the scale of prey diversity is much smaller (range = 0.5-0.75), and thus an increase in of 2.6 times per 1 unit prey diversity is not a meaningful increase.
We also see that we get a warning “2 observations deleted due to missingness”, due to two of our available locations being found on land, where we do not have prey diversity values. This could be mitigated by, for example, generating more than 10 (i.e., 15) available locations per used location, omitting the samples without prey diversity values, and then randomly selecting 10 available locations per used location of those that have values for prey diversity.
Now that we have maximum-likelihood estimates, we will demonstrate how the estimated coefficient for step length can be used as a modifier of the step length shape and scale of the gamma distributions (Avgar et al. 2015). See Avgar et al. (2015) for details.
# grab observation step length parameter
obs_sl_par <- fit_distr(data_obs_steps$sl_, "gamma") # by default, random_steps() uses fit_distr(x$sl_, "gamma")
# grab observed data step length shape
b1 <- obs_sl_par$params$shape
# grab observed data step length scale
b2 <- obs_sl_par$params$scale
# shape after prey diversity selection is included
ssf1_gamma_shape <- b1 + summary(ssf1)$coef["log_sl",1]
# notice that this modified shape is a bit different from the observed shape
ssf1_gamma_shape # modified shape## [1] 1.39735b1 # observed shape, a bit different## [1] 1.355159# we can also modify scale to include prey diversity
ssf1_gamma_scale <- 1/(1/b2 - summary(ssf1)$coef["sl_",1])
# notice that it is a bit different from the observed scale
ssf1_gamma_scale # modified scale## [1] 8452.295b2 # observed scale, a bit different## [1] 8916.152# Since by definition, the mean of the gamma is shape*scale, we will show the different gamma means:
# Mean step length from gamma parametric with observed step length
b1*b2## [1] 12082.8# Mean step length from gamma modified by prey diversity.
ssf1_gamma_shape*ssf1_gamma_scale## [1] 11810.82Next we will interpret ssf2.
# see model summary
summary(ssf2)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ sl_ + preydiv *
## log_sl + cos_ta + strata(step_id_), data = data, model = TRUE,
## method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sl_ -1.068e-05 1.000e+00 1.657e-05 -0.645 0.519
## preydiv 3.150e+01 4.770e+13 2.309e+01 1.364 0.173
## log_sl 1.981e+00 7.250e+00 1.420e+00 1.395 0.163
## cos_ta 2.138e-02 1.022e+00 1.311e-01 0.163 0.871
## preydiv:log_sl -3.123e+00 4.401e-02 2.253e+00 -1.386 0.166
##
## exp(coef) exp(-coef) lower .95 upper .95
## sl_ 1.000e+00 1.000e+00 1.000e+00 1.000e+00
## preydiv 4.770e+13 2.097e-14 1.057e-06 2.152e+33
## log_sl 7.250e+00 1.379e-01 4.483e-01 1.172e+02
## cos_ta 1.022e+00 9.788e-01 7.900e-01 1.321e+00
## preydiv:log_sl 4.401e-02 2.272e+01 5.318e-04 3.643e+00
##
## Concordance= 0.554 (se = 0.029 )
## Likelihood ratio test= 2.24 on 5 df, p=0.8
## Wald test = 2.02 on 5 df, p=0.8
## Score (logrank) test = 2.12 on 5 df, p=0.8Similar to ssf1, we don’t see any significant relationships. Note that we also included a term for preydiv*log_sl. Since we have this interaction, preydiv and log_sl cannot be interpreted independently, thus we interpret the interaction. The interaction has a negative coefficient, suggesting that step lengths decreases with increased prey diversity, which makes sense intuitively, although note that the interaction is not significant.
3.3 HMM
3.3.1 2-state HMMs
Finally, we will fit the HMMs using momentuHMM (McClintock and Michelot, 2018). In preparation, we define initial parameters, and then update the parameters using our fit model, to ultimately fit a more refined model. We will first fit a two-state model.
This HMM example is parameterized to allow the covariate, prey diversity, influence the state transition probability. We can also allow covariates to influence the observation probabilities; these covariates are typically abiotic, but we show a brief example of how one would fit such a model below our main HMM.
## prepare the data
data_hmm <- seal %>%
make_track(lon, lat, date) %>%
extract_covariates(fish_raster) %>%
mutate(ID = 1,
x = x_,
y = y_,
date = t_) %>%
dplyr::select(ID, x, y, date, preydiv) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y")) %>%
st_set_crs("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") %>%
st_transform("+proj=longlat +datum=WGS84") %>%
mutate(long = unlist(map(geometry,1)),
lat = unlist(map(geometry,2))) %>%
st_drop_geometry()
# calculate step length
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("long", "lat"), type = "LL", covNames = c("preydiv"))3.3.1.1 Basic 2-state HMM
Define the starting parameters.
# define parameters
nbStates <- 2 # number of states
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
dist <- list(step=stepDist,angle=angleDist)
mu0 <- c(8, 38) # mean step length for each state
sigma0 <- c(8, 8) # sd step length for each state
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.35) # turning angle for each stateFit the two-state HMM.
set.seed(2023)
hmm_trans_2 <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
stateNames = c("Slow movement", "Moderate movement"),
dist=dist,
Par0=list(step=stepPar0,angle=kappa0),
retryFits = 2)## Attempting to improve fit using random perturbation. Press 'esc' to force exit from 'fitHMM'
##
Attempt 1 of 2 -- current log-likelihood value: -702.4317
Attempt 2 of 2 -- current log-likelihood value: -702.4317As starting values can influence the model fit, it is advisable to fit the model with multiple starting values. This is what the argument retryFits achieves. To limit computation time of the tutorial, it is limited to 2, but a larger number could be used in a real analysis.
3.3.1.2 2-states with covariate effects on transition probabilities
Here, we will explore adding effect of covariates on transition probabilities. Define transition probability model.
formula = ~ preydiv # identify covariatesRetrieve parameters to refine the model.
Par0_hmm_trans_2 <- momentuHMM::getPar0(hmm_trans_2, formula=formula)Note: when new covariates are added to using getPar0(), the intercept coefficients are taken from the previous model, while initial value for slope coefficients are set to 0. We can see this by comparing the estimated transition probability parameters from the previous model to the initial parameters prepared for the next model with covariates.
# transition probability matrix (TPM) parameters from hmm1_km
hmm_trans_2$mle$beta## 1 -> 2 2 -> 1
## (Intercept) -3.312105 -0.5739618# initial parameters for hmm2_km with effect of preydiv on TPM
Par0_hmm_trans_2$beta## 1 -> 2 2 -> 1
## (Intercept) -3.312105 -0.5739618
## preydiv 0.000000 0.0000000Fit a refined HMM with effect of prey diversity on transition probability using estimated parameters from the initial two-state HMM.
set.seed(2023)
hmm_trans_2_refined <- momentuHMM::fitHMM(data=data_hmm, nbStates=2,
stateNames = c("Slow movement", "Moderate movement"),
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm_trans_2$Par,
delta0 = Par0_hmm_trans_2$delta,
beta0 = Par0_hmm_trans_2$beta,
formula=formula, retryFits = 2)## Attempting to improve fit using random perturbation. Press 'esc' to force exit from 'fitHMM'
##
Attempt 1 of 2 -- current log-likelihood value: -701.7688
Attempt 2 of 2 -- current log-likelihood value: -701.7688Visualize the results.
plot(hmm_trans_2_refined)## Decoding state sequence... DONE
We can see from the step-length histogram and the map that it doesn’t look like this HMM is capturing the short and moderate step length movements as their own states. We can look at the pseudo-residuals of the model to confirm whether this is the case.
plotPR(hmm_trans_2_refined)## Computing pseudo-residuals...
We can see in the Q-Q plot for step length (top row, middle plot), that the observed step lengths (red points) are have a higher sample quantile than predicted by the model (shaded grey area) around the theoretical quantile around 1.5. This can be interpreted as the model predicting more steps around the 1.5th theoretical quantile than there are, or in other words, that the speed of steps around the 1.5th quantile are faster than the model predicts. There is also some autocorrelation in the step length (top row, right plot) above the dotted blue line, suggesting that the latent states do not explain all the variation in observed step length. The residuals may be due to a missing state or missing effect of covariates on observation process. We will explore both of these.
3.3.2 3-state HMMs
Next, we will explore adding a third state to the HMM. We start with a basic HMM without any covariates, then explore adding covariates to transition and observation probabilities.
3.3.2.1 Basic 3-state HMM
First, we define starting parameters for a three-state HMM.
# define parameters
nbStates <- 3 # number of states
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
mu0 <- c(5, 12, 38) # mean step length for each state
sigma0 <- c(3, 5, 8) # sd step length for each state
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5) # turning angle for each stateFit the three-state HMM.
set.seed(2023)
hmm_3 <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
stateNames = c("Slow movement", "Moderate movement", "Fast movement"),
dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar0,angle=kappa0), retryFits = 2)## Attempting to improve fit using random perturbation. Press 'esc' to force exit from 'fitHMM'
##
Attempt 1 of 2 -- current log-likelihood value: -679.384
Attempt 2 of 2 -- current log-likelihood value: -679.384Explore model.
plot(hmm_3, ask = F)## Decoding state sequence... DONE
plotPR(hmm_3)## Computing pseudo-residuals...
The predicted state map seems to adequately capture slow, moderate, and fast movement. Importantly, the ACF plots are significantly improved from previous 2-state models. Therefore, the remainder of our analysis will be built on the 3-state model.
3.3.2.2 3-states with covariate effects on transition probabilities
Here, we will explore adding effects of covariates on the transition probabilities. Note that this is the most typical application of covariates to HMM, and this model is the main model we focus on in the main paper.
First, define transition probability model, and retrieve parameters to refine the model.
formula = ~ preydiv # identify covariates
Par0_hmm3 <- momentuHMM::getPar0(hmm_3, formula=formula)Fit a refined HMM with the parameters from the initial three-state HMM.
set.seed(2023)
hmm_trans_3 <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
stateNames = c("Slow movement", "Moderate movement", "Fast movement"),
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm3$Par,
delta0 = Par0_hmm3$delta,
beta0 = Par0_hmm3$beta,
formula=formula, retryFits = 2)## Attempting to improve fit using random perturbation. Press 'esc' to force exit from 'fitHMM'
##
Attempt 1 of 2 -- current log-likelihood value: -675.4028
Attempt 2 of 2 -- current log-likelihood value: -675.4028We can view the parameters and regression coefficients for the transition probabilities.
hmm_trans_3$CIbeta## $step
## $step$est
## mean_1:(Intercept) mean_2:(Intercept) mean_3:(Intercept) sd_1:(Intercept)
## [1,] 1.544323 2.669887 3.68362 1.144003
## sd_2:(Intercept) sd_3:(Intercept)
## [1,] 1.580641 1.424134
##
## $step$se
## mean_1:(Intercept) mean_2:(Intercept) mean_3:(Intercept) sd_1:(Intercept)
## [1,] 0.141362 0.05561547 0.02967768 0.2033532
## sd_2:(Intercept) sd_3:(Intercept)
## [1,] 0.132132 0.2030604
##
## $step$lower
## mean_1:(Intercept) mean_2:(Intercept) mean_3:(Intercept) sd_1:(Intercept)
## [1,] 1.267259 2.560883 3.625453 0.7454379
## sd_2:(Intercept) sd_3:(Intercept)
## [1,] 1.321667 1.026143
##
## $step$upper
## mean_1:(Intercept) mean_2:(Intercept) mean_3:(Intercept) sd_1:(Intercept)
## [1,] 1.821388 2.778892 3.741787 1.542568
## sd_2:(Intercept) sd_3:(Intercept)
## [1,] 1.839615 1.822125
##
##
## $angle
## $angle$est
## concentration_1:(Intercept) concentration_2:(Intercept)
## [1,] -1.132571 -0.3513391
## concentration_3:(Intercept)
## [1,] 0.3820636
##
## $angle$se
## concentration_1:(Intercept) concentration_2:(Intercept)
## [1,] 0.5790591 0.3212703
## concentration_3:(Intercept)
## [1,] 0.37962
##
## $angle$lower
## concentration_1:(Intercept) concentration_2:(Intercept)
## [1,] -2.267506 -0.9810172
## concentration_3:(Intercept)
## [1,] -0.3619779
##
## $angle$upper
## concentration_1:(Intercept) concentration_2:(Intercept)
## [1,] 0.002364172 0.2783391
## concentration_3:(Intercept)
## [1,] 1.126105
##
##
## $beta
## $beta$est
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) 11.95759 -28.98618 0.8889846 -6.008464 52.18805 -12.38337
## preydiv -22.38610 0.00000 -3.8455964 6.298986 -97.11339 19.08809
##
## $beta$se
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) 6.490798 1.901943e-05 6.070081 7.694364 0.8924444 9.787294
## preydiv 10.919036 NaN 9.973558 12.461283 0.5617634 15.703774
##
## $beta$lower
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -0.7641397 -28.98622 -11.00815 -21.08914 50.43889 -31.56611
## preydiv -43.7870196 NaN -23.39341 -18.12468 -98.21442 -11.69075
##
## $beta$upper
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) 24.679321 -28.98614 12.78612 9.072212 53.93721 6.799378
## preydiv -0.985184 NaN 15.70222 30.722652 -96.01235 49.866916
##
##
## $delta
## $delta$est
## Moderate movement Fast movement
## (Intercept) 13.03817 -3.673484
##
## $delta$se
## Moderate movement Fast movement
## (Intercept) 0.002596365 0.01153526
##
## $delta$lower
## Moderate movement Fast movement
## (Intercept) 13.03308 -3.696092
##
## $delta$upper
## Moderate movement Fast movement
## (Intercept) 13.04326 -3.650875It’s helpful to visualize the results for easier interpretation of the parameters and regression coefficients.
plot(hmm_trans_3)## Decoding state sequence... DONE
plotPR(hmm_trans_3)## Computing pseudo-residuals...
We can see that this three-state HMM does better than the two-state HMM. We can see that those short and moderate step-lengths are captured in state 1 and 2, respectively, and the map of the decoded states appears to group different looking movement patterns as different states. Further, the pseudo-residual plot for step length look pretty good, with almost all of the observed points in the Q-Q plot falling within the predicted Q-Q bands.
There is no apparent effect of prey diversity on transitions from slow movement or moderate movement. However, as prey diversity increases, there appears to be a decrease in the transition probability from travel to slow movement, an increase in transition probability from fast to moderate movement, and a peak in remaining within the travelling state when prey diversity is around 0.57.
statianary_est <- plotStationary(hmm_trans_3, plotCI = TRUE, return = T) # returns a list of data frames with estimated stationary state distribution and confidence intervals (useful for plotting with ggplot)
# combine into long data.frame as required by ggplot
statianary_est <- dplyr::bind_rows(statianary_est[[1]]) %>%
mutate(state = rep(hmm_trans_3$stateNames, each = 101))
# set colors
colours.states <- c("#99DDB6", "#539D9C", "#312C66")
# generate plot
ggplot(statianary_est, aes(x = cov, y = est, fill = state)) +
geom_ribbon(aes(ymin = est-se, ymax = est+se), alpha = 0.2) +
geom_line(aes(col = state)) +
scale_colour_manual(values = colours.states) +
scale_fill_manual(values = colours.states) +
ylab("State Probability") +
xlab("Prey Diversity") +
theme_minimal()
This plot shows that the probability of being in the slow behavior, which is often of importance for conservation because it is thought to be associated with foraging, had a positive relationship with prey diversity, the probability of being in the moderate speed behavior decreases with prey diversity, and the probability of being in fast behavior remained relatively constant with prey diversity, with a potential peak in low-medium prey diversity.
3.3.2.3 3-states with covariate effects on observation probabilities
Here we allow prey diversity to influence the step length and turning angle in the observation probabilities. The covariates allowed to influence the observation probabilities are typically abiotic (e.g., snow depth), but for the sake of this tutorial, we will fit this model using prey diversity.
First we define the formulas. Here, we will focus on the effects of prey diversity on step length mean and turning angle concentration.
distFormula <- ~ preydiv
stepDM <- list(mean = distFormula, sd = ~ 1)
angleDM <- list(concentration = distFormula)
DM <- list(step = stepDM, angle = angleDM)Next we grab the parameters from our refined three-state HMM and update the model to include a design matrix (DM) for effects of covariates on distributions.
Par0_hmm_obs_3 <- getPar0(hmm_trans_3, DM = DM)Just as when we added covariates to the transition probability, when adding new covariates to observation probabilities using getPar0(), the intercept coefficients are taken from the previous model and initial slope coefficients are set to 0. We can see this by comparing the previously estimated parameters to the initial parameters for the next model with covariates. Note that when DM is not included in the model, momentuHMM assumes initial parameters are on the natural scale (i.e., on the scale of how observations are measured). However, when DM is included, then initial parameters must be provided on the working scale (i.e., the scale in which the optimization occurs, between negative and positive infinity). To convert working-scale parameters to the natural scale, we must apply the inverse of the link function associated with that parameter (Table 2 in the vignette for momentuHMM). In the case of mean and standard deviation associated with the gamma distribution used for step length, the link function is log(), therefore we must apply exp() to convert the working scale parameters returned by getPar0() to the natural scale.
Fit the model.
hmm_obs_3 <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
stateNames = c("Slow movement", "Moderate movement", "Fast movement"),
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm_obs_3$Par,
DM = DM, retryFits = 2)## =======================================================================## Fitting HMM with 3 states and 2 data streams## -----------------------------------------------------------------------## step ~ gamma(mean=~preydiv, sd=~1)## angle ~ vm(concentration=~preydiv)##
## Transition probability matrix formula: ~1##
## Initial distribution formula: ~1## =======================================================================## Attempting to improve fit using random perturbation. Press 'esc' to force exit from 'fitHMM'
##
Attempt 1 of 2 -- current log-likelihood value: -675.2531
Attempt 2 of 2 -- current log-likelihood value: -675.2531##
## DONEVisualize the results.
plot(hmm_obs_3, ask = F, plotCI = T)## Decoding state sequence... DONE
The plots show that as prey diversity increases, the step length mean appears to decrease during slow movement. For the moderate movement, step length appears to increase as prey diversity increases. Likely due to low frequency of travelling (fast movement) state, the confidence intervals of effect of prey diversity on step length mean are unreliable.
There is no significant effect of prey diversity in turn angle concentration on any state. While not significant in our data set, turn angle concentration may decrease as prey diversity increases during slow movement (i.e., becomes more tortuous), and may be worth exploring with larger sample sizes.
View parameters.
hmm_obs_3$CIbeta## $step
## $step$est
## mean_1:(Intercept) mean_1:preydiv mean_2:(Intercept) mean_2:preydiv
## [1,] 2.690414 -1.957697 1.115528 2.495038
## mean_3:(Intercept) mean_3:preydiv sd_1:(Intercept) sd_2:(Intercept)
## [1,] 3.815786 -0.2290479 1.082236 1.610697
## sd_3:(Intercept)
## [1,] 1.454613
##
## $step$se
## mean_1:(Intercept) mean_1:preydiv mean_2:(Intercept) mean_2:preydiv
## [1,] 1.317648 2.125084 0.6738791 1.09184
## mean_3:(Intercept) mean_3:preydiv sd_1:(Intercept) sd_2:(Intercept)
## [1,] 0.4106169 0.6782186 0.2022278 0.1188858
## sd_3:(Intercept)
## [1,] 0.2406382
##
## $step$lower
## mean_1:(Intercept) mean_1:preydiv mean_2:(Intercept) mean_2:preydiv
## [1,] 0.1078717 -6.122784 -0.2052512 0.3550701
## mean_3:(Intercept) mean_3:preydiv sd_1:(Intercept) sd_2:(Intercept)
## [1,] 3.010991 -1.558332 0.6858769 1.377685
## sd_3:(Intercept)
## [1,] 0.9829711
##
## $step$upper
## mean_1:(Intercept) mean_1:preydiv mean_2:(Intercept) mean_2:preydiv
## [1,] 5.272957 2.207391 2.436306 4.635005
## mean_3:(Intercept) mean_3:preydiv sd_1:(Intercept) sd_2:(Intercept)
## [1,] 4.62058 1.100236 1.478595 1.843709
## sd_3:(Intercept)
## [1,] 1.926256
##
##
## $angle
## $angle$est
## concentration_1:(Intercept) concentration_1:preydiv
## [1,] 8.037653 -15.01587
## concentration_2:(Intercept) concentration_2:preydiv
## [1,] -3.28635 4.793973
## concentration_3:(Intercept) concentration_3:preydiv
## [1,] 3.344218 -4.922257
##
## $angle$se
## concentration_1:(Intercept) concentration_1:preydiv
## [1,] 7.428896 12.71922
## concentration_2:(Intercept) concentration_2:preydiv
## [1,] 4.843114 7.895969
## concentration_3:(Intercept) concentration_3:preydiv
## [1,] 4.099148 6.96634
##
## $angle$lower
## concentration_1:(Intercept) concentration_1:preydiv
## [1,] -6.522717 -39.94508
## concentration_2:(Intercept) concentration_2:preydiv
## [1,] -12.77868 -10.68184
## concentration_3:(Intercept) concentration_3:preydiv
## [1,] -4.689965 -18.57603
##
## $angle$upper
## concentration_1:(Intercept) concentration_1:preydiv
## [1,] 22.59802 9.913337
## concentration_2:(Intercept) concentration_2:preydiv
## [1,] 6.205979 20.26979
## concentration_3:(Intercept) concentration_3:preydiv
## [1,] 11.3784 8.731517
##
##
## $beta
## $beta$est
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -1.709311 -18.73859 -1.691819 -2.287285 -16.34934 -0.4409921
##
## $beta$se
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) 0.4740712 NaN 0.4488078 0.5062071 NaN 0.5859873
##
## $beta$lower
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -2.638473 NaN -2.571467 -3.279433 NaN -1.589506
##
## $beta$upper
## 1 -> 2 1 -> 3 2 -> 1 2 -> 3 3 -> 1 3 -> 2
## (Intercept) -0.7801483 NaN -0.8121724 -1.295137 NaN 0.7075218
##
##
## $delta
## $delta$est
## Moderate movement Fast movement
## (Intercept) 16.82743 -4.223506
##
## $delta$se
## Moderate movement Fast movement
## (Intercept) NaN NaN
##
## $delta$lower
## Moderate movement Fast movement
## (Intercept) NaN NaN
##
## $delta$upper
## Moderate movement Fast movement
## (Intercept) NaN NaNWe see that when including covariates affecting the observation distributions in a 3-state model, step length mean decreases as prey diversity increases for both the slow movement and fast movement states, where as the step length mean increase as prey diversity increases for the moderate movement state. Likewise, turning angle concentration decreases as prey diversity increases for both the slow movement and fast movement states, where as the turning angle concentration increase as prey diversity increases for the moderate movement state. In other words, as prey diversity increases, behaviors become slower (or faster) and more tortuous.
4 Plots of predicted relationships with the covariate
Here we are plotting a model’s estimated relationship between the resource covariate and probability of selection can be useful for general ecological inference. We will calculate the log of the relative selection strength (log-RSS) for each selection function model. The log-RSS is a measure of how likely a location (for the RSF) or step (for SSFs) is to end in a proposed location (x1) to a single reference location (x2, the mean prey diversity), where zero indicates no preference, >1 indicates selection, and <1 indicates avoidance (Avgar et al. 2017, Fieberg et al. 2021).
First, we prepare a dataframe to predict on.
# prep the fish data
newfish <- fish_raster %>%
terra::as.data.frame(xy = TRUE) %>%
filter(x > 100000 & x < 600000 & y > -550000 & y < 0)4.1 RSF
Since the RSF does not incorporate movement, we will calculate the log-RSS of the movement-free habitat selection kernel. This is easily done using log_rss() from amt.
First, we make a base dataframe to create x1 and x2 from. The values of log_sl and cos_ta do not matter, but we need populated columns in order for the log_rss function to work.
base <- newfish %>%
mutate(sl_ = mean(data_ssf$sl_),
log_sl = mean(data_ssf$log_sl),
cos_ta = cos(1))
# x1 is our base dataframe
x1 <- base Next, we modify the base data frame, where prey diversity is held at its mean.
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))Now we will apply log_rss() to each row. Since log_rss() only assessed one location relative to a reference point, we will use lapply to iterate through all locations.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(rsf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# combine rows
res1 <- dplyr::bind_rows(log_rss_list)Visualize results.
# plot
line_rsf <- ggplot(res1, aes(x = preydiv_x1, y = (log_rss))) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_line(linewidth = 1, color = "#F8D59F",linetype = 2) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.4, fill = "#F8D59F") +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()
line_rsf
We see a positive relationship between prey diversity and log-RSS, which suggests that our seal is more likely to be present in areas with higher prey diversity than areas with low prey diversity. Note the bowtie shape of the standard error is due to log-RSS calculating selection strength relative to a starting step, in this case, where prey diversity is set to its mean.
4.2 SSF
We can also calculate the log-RSS for our SSFs following the same workflow as for the RSF. Since log_rss() passes to predict(), it is important that the fit SSF includes model = TRUE.
First, we will calculate log-RSS for ssf1. Again, since log_rss() only assesses one location relative to a reference point, we will use lapply to iterate through all locations, and bind the rows together after so that we have a single data frame for plotting.
## log-RSS prediction for ssf1
# apply log_rss() to each row
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# combine rows
res2 <- dplyr::bind_rows(log_rss_list) %>%
mutate(Speed = "without int.")Visualize results.
line_ssf1 <- ggplot() +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
geom_line(data = res2, aes(x = preydiv_x1, y = log_rss, color = Speed, group = Speed), size = 1, linetype = 3, color = "black") +
geom_ribbon(data = res2, aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = Speed, group = Speed), alpha = 0.3, fill = "black") +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.line_ssf1
We see no relationship between prey diversity and log-RSS, which suggests that our seal is similarily likely to be present in areas with higher prey diversity and areas with low prey diversity.
Now we will calculate log-RSS for ssf2. We will estimate the log-RSS for three different step-lengths (slow, moderate, fast). We set these speeds as the 25th, 50th, and 75th percentiles of step-length, then loop the log-RSS for each speed. First, identify what the percentiles are.
# determine the 25th (slow), 50th (moderate), and 75th (fast) percentiles of step-length
nums <- seal %>%
make_track(lon, lat, date) %>%
steps %>%
mutate(log_sl = log(sl_)) %>%
dplyr::reframe(quants = quantile(log_sl, c(0.25, 0.5, 0.75))) %>%
pull()Now apply log-RSS to each row, for each speed (step length percentile).
# set-up to run function for each speed
results_ssf2 <- lapply(nums, function(j) {
x1$log_sl <- j
# calculate log-RSS
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf2, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# bind rows within each speed's prediction
res3 <- dplyr::bind_rows(log_rss_list)
} )Bind the output together, and name the step lengths based on how long they are, so that we can interpret them in the subsequent plot as different speeds.
results_ssf2 <- dplyr::bind_rows(results_ssf2) %>%
mutate(log_sl_x1 = as.factor(round(log_sl_x1,1)),
Speed = dplyr::case_when(as.factor(log_sl_x1) == '8.4' ~ "Slow",
as.factor(log_sl_x1) == "9.2" ~ "Moderate",
as.factor(log_sl_x1) == "9.6" ~ "Fast"))Visualize results.
line_ssf2 <- ggplot(results_ssf2, aes(x = preydiv_x1, y = (log_rss))) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_line(size = 1, aes(color = Speed, group = Speed, linetype = Speed)) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = Speed, group = Speed), alpha = 0.3) +
scale_colour_manual(values=c(colours.states,"#000000")) +
scale_fill_manual(values=c(colours.states,"#000000")) +
scale_linetype_manual(values = c("solid", "solid", "solid")) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()
line_ssf2
We can see a weak positive relationship between prey diversity and log-RSS, however, with confidence intervals (shading) covering zero, suggesting no significant relationship between prey diversity and log-RSS. We also see that the confidence intervals cover each other, suggesting that different speeds do not have different relationships between prey diversity and log-RSS.
Typically these models would be interpreted independently. But it’s worth noting that while the effects are minimal and the confidence intervals overlap, when comparing ssf1 and ssf2, ssf2 provides more information about the animal’s relationship with prey diversity.
ggplot(results_ssf2, aes(x = preydiv_x1, y = (log_rss))) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_line(size = 1, aes(color = Speed, group = Speed, linetype = Speed)) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = Speed, group = Speed), alpha = 0.3) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal() +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_line(data = res2, aes(x = preydiv_x1, y = log_rss, color = Speed, group = Speed, linetype = Speed), size = 1) +
geom_ribbon(data = res2, aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = Speed, group = Speed), alpha = 0.2) +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66", "black"),
labels = c("Fast", "Moderate", "Slow", "Not consid."),
name = "Speed") +
scale_fill_manual(values = c("#99DDB6", "#539D9C", "#312C66", "black"),
labels = c("Fast", "Moderate", "Slow", "Not consid."),
name = "Speed") +
scale_linetype_manual(values = c("solid", "solid", "solid", "dotted"),
labels = c("Fast", "Moderate", "Slow", "Not consid."),
name = "Speed") 
Here, we see that adding the interaction between prey diversity and step length did not better explain our data or provide additional ecological insight about our seal.
4.3 HMM
Exploring the insight on the relationship between our seal and prey diversity is fundamentally different for the HMM. Here, we will grab and plot the stationary state probabilities. This is easily done using plotStationary() from momentuHMM.
# grab stationary probabilities
ps <- momentuHMM::plotStationary(hmm_trans_3, plotCI= TRUE, return = TRUE)# grab values for data frame
state1 <- ps$preydiv$'Slow movement' %>% mutate(state = 1)
state2 <- ps$preydiv$'Moderate movement' %>% mutate(state = 2)
state3 <- ps$preydiv$'Fast movement' %>% mutate(state = 3)
# bind to one data frame
pdat <- rbind(state1, state2, state3) %>%
mutate(state = as.character(state))Visualize results.
line_hmm <- ggplot() +
geom_line(data = pdat, aes(x = cov, y = est, color = state)) +
geom_ribbon(data = pdat, aes(x=cov, y=est, ymax=est+se, ymin=est-se, fill = state),
alpha = 0.4, show.legend = TRUE) +
ylab("State probabilty") +
xlab("Prey diversity") +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Slow movement", "Moderate movement", "Fast movement")) +
scale_fill_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Slow movement", "Moderate movement", "Fast movement")) +
theme_minimal()
line_hmm
Here we can see that each state has a different relationship between the stationary state probability and prey diversity. The slow movement state has a positive relationship between stationary state probability and prey diversity, the moderate movement state has a negative relationship with prey diversity, and the fast movement state does not appear to have a directional relationship with prey diversity.
5 Prediction maps
Now we will estimate the utilization distributions from each model to demonstrate how differences in the relationships with a covariate can results in vastly different spatial patterns. The utilization distribution is defined as the two-dimensional relative frequency distribution of space use of an animal (Van Winkle 1975). This is a simple calculation for the RSF, where we multiply the model coefficient with the resource (prey diversity), exponentiate (since it is a logistic regression), and normalize the estimate. The calculations are more complex for the SSFs since they are conditional models that integrate the movement process. Thus, for the SSFs we calculate the steady-state utilization distribution (SSUD), which is the long-term expectation of the space-use distribution across the landscape (Signer et al. 2017). amt has functions to estimate the SSUD.
5.1 RSF
We can predict the estimated probability of use from the RSF by hand. First, grab the model coefficients and predict for each cell.
# grab model coefficients
modcoef <- summary(rsf1)$coef
# prediction for each cell
x <- exp(modcoef[2] * newfish$preydiv)We will normalize the results next.
# range fn
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
# set the range from zero to one
newfish$rsf_prediction <- range01(x)Visualize results.
map_rsf <- ggplot() +
geom_tile(data = newfish, aes(x = x,y = y, fill = rsf_prediction)) +
scale_fill_viridis(option = "mako", name = "RSF prediction",limits = c(0,1)) +
geom_sf(data = nat_trans, fill = "grey80", color = "white") +
coord_sf(xlim = ext(fish_raster)[1:2], ylim = ext(fish_raster)[3:4], expand = F) +
theme_void()
map_rsf
5.2 SSF
We can use the amt functions to estimate the SSUDs from the simple SSF that does not allow prey diversity to affect the movement kernel. When estimating the SSUDs, the functions require that covariates are specifically labelled as occuring at the end of the step (i.e., preydiv_end). This is easily implemented by specifying where = "both" when using extract_covariates, which will then produce a column for preydiv_start and preydiv_end. By default, as in our initial fitting, extract_covariates extracts data for the end of the step, but does not record the _start or _end suffix. Thus, we are re-extracting covariates for ease throughout the generation of the SSUDs.
# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date) %>%
steps() %>%
random_steps() %>%
arrange(case_) %>%
amt::extract_covariates(fish_raster, where = "both")
# fit SSF1 model
m1 <- data_ssf |>
fit_clogit(case_ ~ preydiv_end + sl_ + log(sl_) + cos(ta_) +
strata(step_id_))Now we will set the starting position for the simulation.
start <- make_start((seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date))[1,])We can also set constants for how many steps we want to simulate. Additionally, we set a burn-in number of steps, where we later remove this number of steps from the simulated path to reduce the influence of the starting point.
n_steps = 1e4 # number of steps
burnin <- n_steps/50 # number of steps to remove for the burn-inGenerate the redistribution kernel.
k1 <- redistribution_kernel(m1, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 1,
n.control = n_steps)Simulate the path. This function will likely take about five minutes on most laptops. As such, I have added the resulting simulated paths /data/p1_1e4 Feel free to uncomment the read.csv line to load it in instead of running the simulate_path if prefered. Additionally, I have also ran the simulation with 1e5 steps. This takes about 12 hours to run, so we will not use those data here, but this can be loaded as /data/p1_1e4. Our plots in the main paper use 1e5 steps.
set.seed(2023)
p1 <- amt::simulate_path(x = k1, n = n_steps, start = start, verbose = TRUE)
#p1 <- read.csv(here("data/p1_1e4.csv")) %>% as_tibble()Remove the burn-in.
p1_burnt <- p1 %>% dplyr::slice(-c(1:burnin))Visualize the simulated path.
ssf_track_1 <- fishmap +
geom_sf(data = nat_trans, fill = "grey80", color = "white") +
geom_point(data = p1_burnt, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p1_burnt, aes(x = x_,y = y_)) +
coord_sf(xlim = ext(fish_raster)[1:2], ylim = ext(fish_raster)[3:4], expand = F) +
theme_void()## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.ssf_track_1
We will use this track to estimate the SSUD. We will convert the estimated SSUD to a data frame for easy plotting.
uds_ssf1 <- tibble(rep = 1:nrow(p1_burnt),
x_ = p1_burnt$x_, y_ = p1_burnt$y_,
t_ = p1_burnt$t_, dt = p1_burnt$dt) %>%
filter(!is.na(x_)) %>%
make_track(x_, y_) %>%
hr_kde(trast = fish_raster, which_min = "global") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)Visualize the results.
map_ssf1 <- ggplot() +
geom_tile(data = uds_ssf1, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF1 prediction") +
geom_sf(data = nat_trans, fill = "grey80", color = "white") +
coord_sf(xlim = ext(fish_raster)[1:2], ylim = ext(fish_raster)[3:4], expand = F) +
theme_void()
map_ssf1
The map shows that probability of use is fairly generally distributed in space.
We will follow the same steps to generate a SSUD for the SSF that allows prey diversity to affect the movement kernel.
m2 <- data_ssf %>%
fit_clogit(case_ ~ preydiv_end * log(sl_) + sl_ + cos(ta_) +
strata(step_id_))Generate the redistribution kernel.
set.seed(2023)
k2 <- redistribution_kernel(m2, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 1,
n.control = n_steps)Simulate the path. Similar to SSF1, this function may take some time to run. I have stored the simulated path with 1e5 locations in data/p2_1e4, or for the full simulation, data/p2_1e5, for loading in here if you prefer not to run it on your own computer.
set.seed(2023)
p2 <- amt::simulate_path(x = k2, n = n_steps, start = start, verbose = TRUE)
#p2 <- read.csv(here("data/p2_1e4.csv")) %>% as_tibble()Remove burn-in.
p2_burnt <- p2 %>% dplyr::slice(-c(1:burnin))Visualize simulated track.
ssf_track_2 <- fishmap +
geom_sf(data = nat_trans, fill = "grey80", color = "white") +
geom_point(data = p2_burnt, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p2_burnt, aes(x = x_,y = y_)) +
coord_sf(xlim = ext(fish_raster)[1:2], ylim = ext(fish_raster)[3:4], expand = F) +
theme_void()## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.ssf_track_2
Estimate the SSUD.
uds_ssf2 <- tibble(rep = 1:nrow(p2_burnt),
x_ = p2_burnt$x_, y_ = p2_burnt$y_,
t_ = p2_burnt$t_, dt = p2_burnt$dt) |>
filter(!is.na(x_)) |>
make_track(x_, y_) |>
hr_kde(trast = fish_raster, which_min = "global") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)Visualize the SSUD.
map_ssf2 <- ggplot() +
geom_tile(data = uds_ssf2, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF2 prediction") +
geom_sf(data = nat_trans, fill = "grey80", color = "white") +
coord_sf(xlim = ext(fish_raster)[1:2], ylim = ext(fish_raster)[3:4], expand = F) +
theme_void()
map_ssf2
5.3 HMM
We will first estimate the stationary state probabilities of each state based on prey diversity. This is easily done using the momentuHMM function stationary(). This function predicts the probability of being in each state given a covariate (prey diversity).
x <- as.data.frame(momentuHMM::stationary(hmm_trans_3,
data.frame(preydiv = newfish$preydiv)))
newfish$hmm_state1 <- x$Slow.movement
newfish$hmm_state2 <- x$Moderate.movement
newfish$hmm_state3 <- x$Fast.movementFormat the data for plotting.
newfish_long <- newfish %>%
tidyr::pivot_longer(cols = hmm_state1:hmm_state3,
names_to = "model", values_to = "prediction") %>%
mutate(dplyr::across(model, factor, levels=
c("hmm_state1", "hmm_state2", "hmm_state3")))## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `dplyr::across(...)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))Visualize the results.
map_hmm <- ggplot() +
geom_tile(data = newfish_long, aes(x = x, y = y, fill = prediction)) +
scale_fill_viridis(option = "mako", limits = c(0,1)) +
labs(fill = 'HMM predicted\nprobability') +
geom_sf(data = nat_trans, fill = "grey80", color = "white") +
coord_sf(xlim = ext(fish_raster)[1:2], ylim = ext(fish_raster)[3:4], expand = F) +
facet_wrap(~model, labeller = as_labeller(c('hmm_state1' = "Slow movement",
'hmm_state2' = "Moderate movement",
'hmm_state3' = "Fast movement"))) +
theme_void()
map_hmm
We can see that the slow movement state, which was positively related to prey diversity, showed a similar restricted spatial pattern to the RSF, whereas the moderate movement state, which had a negative relationship with prey diversity, showed the opposite pattern to the RSF prediction map. The fast movement state did not show as much spatial variation, due to its non significant relationship with prey diversity.
6 References
7 Session
sessioninfo::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.1.0 (2021-05-18)
## os macOS Big Sur 10.16
## system x86_64, darwin17.0
## ui X11
## language (EN)
## collate en_CA.UTF-8
## ctype en_CA.UTF-8
## tz America/Vancouver
## date 2024-01-24
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.1.0)
## amt * 0.2.2.0 2023-07-07 [1] Github (jmsigner/amt@3622dbd)
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.1.0)
## backports 1.2.1 2020-12-09 [1] CRAN (R 4.1.0)
## bitops 1.0-7 2021-04-24 [1] CRAN (R 4.1.0)
## boot 1.3-28 2021-05-03 [1] CRAN (R 4.1.0)
## Brobdingnag 1.2-7 2022-02-03 [1] CRAN (R 4.1.2)
## broom 0.7.9 2021-07-27 [1] CRAN (R 4.1.0)
## bslib 0.3.1 2021-10-06 [1] CRAN (R 4.1.0)
## car 3.0-11 2021-06-27 [1] CRAN (R 4.1.0)
## carData 3.0-4 2020-05-22 [1] CRAN (R 4.1.0)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.1.0)
## checkmate 2.0.0 2020-02-06 [1] CRAN (R 4.1.0)
## CircStats 0.2-6 2018-07-01 [1] CRAN (R 4.1.0)
## circular 0.4-94 2022-03-07 [1] CRAN (R 4.1.2)
## class 7.3-19 2021-05-03 [1] CRAN (R 4.1.0)
## classInt 0.4-3 2020-04-07 [1] CRAN (R 4.1.0)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.1.0)
## codetools 0.2-18 2020-11-04 [1] CRAN (R 4.1.0)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.1.0)
## crawl 2.2.1 2018-09-14 [1] CRAN (R 4.1.0)
## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.1.2)
## curl 5.0.0 2023-01-12 [1] CRAN (R 4.1.2)
## data.table 1.14.0 2021-02-21 [1] CRAN (R 4.1.0)
## DBI 1.1.1 2021-01-15 [1] CRAN (R 4.1.0)
## dbplyr 2.1.1 2021-04-06 [1] CRAN (R 4.1.0)
## digest 0.6.32 2023-06-26 [1] CRAN (R 4.1.0)
## doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.1.2)
## doRNG 1.8.2 2020-01-27 [1] CRAN (R 4.1.0)
## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.1.2)
## e1071 1.7-7 2021-05-23 [1] CRAN (R 4.1.0)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.1.0)
## evaluate 0.19 2022-12-13 [1] CRAN (R 4.1.2)
## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.1.0)
## farver 2.1.1 2022-07-06 [1] CRAN (R 4.1.2)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.1.0)
## fitdistrplus 1.1-8 2022-03-10 [1] CRAN (R 4.1.2)
## forcats * 0.5.1 2021-01-27 [1] CRAN (R 4.1.0)
## foreach 1.5.1 2020-10-15 [1] CRAN (R 4.1.0)
## foreign 0.8-81 2020-12-22 [1] CRAN (R 4.1.0)
## fs 1.6.2 2023-04-25 [1] CRAN (R 4.1.0)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.1.2)
## geosphere * 1.5-10 2019-05-26 [1] CRAN (R 4.1.0)
## ggmap 3.0.0 2019-02-05 [1] CRAN (R 4.1.0)
## ggplot2 * 3.4.0 2022-11-04 [1] CRAN (R 4.1.2)
## glue 1.6.2 2022-02-24 [1] CRAN (R 4.1.2)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.1.0)
## gtable 0.3.1 2022-09-01 [1] CRAN (R 4.1.2)
## haven 2.4.3 2021-08-04 [1] CRAN (R 4.1.0)
## here * 1.0.1 2020-12-13 [1] CRAN (R 4.1.0)
## highr 0.9 2021-04-16 [1] CRAN (R 4.1.0)
## hms 1.1.0 2021-05-17 [1] CRAN (R 4.1.0)
## htmltools 0.5.2 2021-08-25 [1] CRAN (R 4.1.0)
## httpuv 1.6.3 2021-09-09 [1] CRAN (R 4.1.0)
## httr 1.4.2 2020-07-20 [1] CRAN (R 4.1.0)
## iterators 1.0.13 2020-10-15 [1] CRAN (R 4.1.0)
## jpeg 0.1-9 2021-07-24 [1] CRAN (R 4.1.0)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.1.0)
## jsonlite 1.8.4 2022-12-06 [1] CRAN (R 4.1.2)
## KernSmooth 2.23-20 2021-05-03 [1] CRAN (R 4.1.0)
## knitr 1.45 2023-10-30 [1] CRAN (R 4.1.0)
## labeling 0.4.2 2020-10-20 [1] CRAN (R 4.1.0)
## later 1.3.0 2021-08-18 [1] CRAN (R 4.1.0)
## lattice 0.20-44 2021-05-02 [1] CRAN (R 4.1.0)
## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.1.2)
## lubridate * 1.9.2 2023-02-10 [1] CRAN (R 4.1.2)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.1.2)
## MASS 7.3-54 2021-05-03 [1] CRAN (R 4.1.0)
## Matrix 1.3-3 2021-05-04 [1] CRAN (R 4.1.0)
## mime 0.11 2021-06-23 [1] CRAN (R 4.1.0)
## modelr 0.1.8 2020-05-19 [1] CRAN (R 4.1.0)
## momentuHMM * 1.5.4 2021-09-03 [1] CRAN (R 4.1.0)
## moveHMM 1.7 2019-05-19 [1] CRAN (R 4.1.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.1.0)
## mvtnorm 1.1-2 2021-06-07 [1] CRAN (R 4.1.0)
## numDeriv 2016.8-1.1 2019-06-06 [1] CRAN (R 4.1.0)
## openxlsx 4.2.4 2021-06-16 [1] CRAN (R 4.1.0)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.1.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.1.0)
## plyr 1.8.7 2022-03-24 [1] CRAN (R 4.1.2)
## png 0.1-7 2013-12-03 [1] CRAN (R 4.1.0)
## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.1.0)
## proxy 0.4-27 2022-06-09 [1] CRAN (R 4.1.0)
## purrr * 1.0.1 2023-01-10 [1] CRAN (R 4.1.2)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.1.0)
## rbibutils 2.2.7 2021-12-07 [1] CRAN (R 4.1.0)
## Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.1.0)
## Rdpack 2.1.4 2022-02-18 [1] CRAN (R 4.1.2)
## readr * 2.1.2 2022-01-30 [1] CRAN (R 4.1.2)
## readxl 1.3.1 2019-03-13 [1] CRAN (R 4.1.0)
## reprex 2.0.2 2022-08-17 [1] CRAN (R 4.1.2)
## rgdal * 1.5-23 2021-02-03 [1] CRAN (R 4.1.0)
## RgoogleMaps 1.4.5.3 2020-02-12 [1] CRAN (R 4.1.0)
## rio 0.5.27 2021-06-21 [1] CRAN (R 4.1.0)
## rjson 0.2.20 2018-06-08 [1] CRAN (R 4.1.0)
## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.1.0)
## rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.1.0)
## rnaturalearth * 0.1.0 2017-03-21 [1] CRAN (R 4.1.0)
## rnaturalearthdata * 0.1.0 2017-02-21 [1] CRAN (R 4.1.0)
## rngtools 1.5.2 2021-09-20 [1] CRAN (R 4.1.0)
## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.1.2)
## rstudioapi 0.13 2020-11-12 [1] CRAN (R 4.1.0)
## rvest 1.0.1 2021-07-26 [1] CRAN (R 4.1.0)
## sass 0.4.0 2021-05-12 [1] CRAN (R 4.1.0)
## scales 1.2.1 2022-08-20 [1] CRAN (R 4.1.2)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 4.1.0)
## sf * 1.0-12 2023-03-19 [1] CRAN (R 4.1.2)
## shiny 1.7.1 2021-10-02 [1] CRAN (R 4.1.0)
## sp * 1.4-5 2021-01-10 [1] CRAN (R 4.1.0)
## stringi 1.7.6 2021-11-29 [1] CRAN (R 4.1.0)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 4.1.0)
## survival 3.2-11 2021-04-26 [1] CRAN (R 4.1.0)
## terra * 1.7-3 2023-01-24 [1] CRAN (R 4.1.2)
## tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.1.2)
## tidyr * 1.2.1 2022-09-08 [1] CRAN (R 4.1.2)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.1.2)
## tidyterra * 0.3.1 2022-11-09 [1] CRAN (R 4.1.2)
## tidyverse * 1.3.1 2021-04-15 [1] CRAN (R 4.1.0)
## timechange 0.2.0 2023-01-11 [1] CRAN (R 4.1.2)
## tzdb 0.1.2 2021-07-20 [1] CRAN (R 4.1.0)
## units 0.7-2 2021-06-08 [1] CRAN (R 4.1.0)
## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.1.0)
## vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.1.0)
## viridis * 0.6.2 2021-10-13 [1] CRAN (R 4.1.0)
## viridisLite * 0.4.2 2023-05-02 [1] CRAN (R 4.1.0)
## withr 2.5.0 2022-03-03 [1] CRAN (R 4.1.2)
## xfun 0.39 2023-04-20 [1] CRAN (R 4.1.0)
## xml2 1.3.2 2020-04-23 [1] CRAN (R 4.1.0)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.1.0)
## yaml 2.2.1 2020-02-01 [1] CRAN (R 4.1.0)
## zip 2.2.0 2021-05-31 [1] CRAN (R 4.1.0)
##
## [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library